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We performed experiments using Aplysia neurons to identify the mechanism underlying the changes 
in the firing patterns in response to temperature changes. When the temperature was gradually 
increased from 11°C to 31°C the firing patterns changed sequentially from the silent state to beating, 
doublets, beating-chaos, bursting-chaos, square-wave bursting, and bursting-oscillation patterns. When 
the temperature was decreased over the same temperature range, these sequential changes in the 
firing patterns reappeared in reverse order. To simulate this entire range of spiking patterns we 
modified nonlinear differential equations that Chay and Lee made using temperature-dependent 
scaling factors. To refine the equations, we also analyzed the spike pattern changes in the presence 
of potassium channel blockers. Based on the solutions of these equations and potassium channel blocker 
experiments, we found that, as temperature increases, the maximum value of the potassium channel 
relaxation time constant, r„(t) increases, but the maximum value of the probabilities of openings for 
activation of the potassium channels, n(t) decreases. Accordingly, the voltage-dependent potassium 
current is likely to play a leading role in the temperature-dependent changes in the firing patterns 
in Aplysia neurons. 
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INTRODUCTION 

It is well known that the general conclusion from other 
studies on the temperature dependence of firing patterns 
of Aplysia neurons was the increase in spike frequency 
[1,2], an increase in resting membrane potential [1,3], and 
a decrease in spike-broadening and spike height [4] with 
increasing temperature. For example, Hakozaki et al. found 
that temperature-induced hyperpolarization was mediated 
in Aplysia neurons by the temperature- sensitive and G-pro- 
tein regulated potassium current as the temperature in- 
creased [5]. AP broadening in Aplysia R15 pacemaker neu- 
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ron resulted from the depolarizing effect of slow volt- 
age-dependent Ca 2+ currents and its peak height was de- 
creased by steady activation of a temperature-sensitive K + 
current [4]. However, the mechanisms of the firing pattern 
changes in response to temperature are still unclear. To 
identify the mechanism of the temperature-induced firing 
pattern change, we examined the effect of temperature on 
the firing pattern by using the neurons of Aplysia and de- 
veloped a computational model to simulate the temper- 
ature-dependent changes in the firing patterns. 

Aplysia neurons have great advantages for electro- 
physiological study. In Aplysia, the nervous system is rela- 
tively simple and the size of individual neuron is so large 
that each neuron can be easily identified based on the loca- 
tion in the ganglia. In addition, Aplysia giant neurons are 
long-lived in vitro and so are apt to stand heat stress for 
a long time. The neurons in the left caudal quarter- ganglion 
(LUQ) of the abdominal ganglion of Aplysia fires spon- 
taneously. We found that many of these neurons in the 
LUQ changed the firing pattern sequentially with temper- 
ature from a silent state to beating, beating-chaos, burst- 
ing-chaos, square-wave bursting, and finally a burst- 

ABBREVIATIONS: TEA, tetraethylammonium chloride; 4-AP, 4-Ami- 
nopyridine; AP, action potential; ISI, the interspike interval; IBI, 
the interburst interval; ASW, artificial seawater; DAQ, data acqui- 
sition; ID, interval duration. 
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ing-oscillation pattern. When the temperature was de- 
creased over the same temperature range, the above seven 
patterns reappeared in reverse order. 

Based on these experimental results, we wanted to com- 
pose equations to simulate all of the seven different spiking 
patterns. However, we found that these various signal pat- 
terns observed between 25°C and 31°C were similar to pic- 
tures published by Chay [6], and Chay and Lee [7]. 
Although in these papers various spiking patterns were de- 
scribed by solving nonlinear differential equations such as 
continuous spiking, doublets, beating-chaos, bursting-chaos, 
square-wave bursting, and slow-wave, the silent phase and 
low frequency beating were not described between 11°C and 
25°C that are observed in our experiments. Furthermore, 
any temperature-dependent scaling factors were not in- 
cluded in these equations, either. Therefore, simply chang- 
ing the variables presented in papers by Chay and Lee [6,7] 
would not be sufficient to simulate these exceptional spik- 
ing patterns. Furthermore it was uncertain whether or not 
the temperature-dependent neurons from the abdominal 
ganglion of A. juliana could be regarded as warm receptors. 
Because cold and warm receptors respond to temperature 
change in one direction only; warm (cold) receptor responds 
to temperature changes during the rising (falling) phase 
but do not respond to temperature changes during the fall- 
ing (rising) phase. Therefore, these neurons are thought to 
be different from warm receptors. However, the effects of 
temperature changes on cold or warm receptors have been 
studied in various species, including cats [8,9], vampire 
bats and mice [10], snakes [11], rabbit [12], and human 
[13]. Although mammalian cold receptors are discussed in 
the Huber-Braun cold receptor model [14,15], it was not 
even able to find a proper biophysical model for warm 
receptors. 

To set up nonlinear differential equations that can simu- 
late all of these various spiking patterns together, it is nec- 
essary to modify the equations derived by Chay and Lee 
[6,7] with temperature-dependent scaling factors; we will 
adopt a temperature-dependent scaling factor for the ionic 
currents and ionic kinetics that are used in the Huber- 
Braun cold receptor model [14,15]. The cell membrane of 
this new model contains sodium, calcium, and potassium 
channels carrying a fast sodium current, In s , fast calcium 
current, Ic*, and a potassium current, Ik, respectively. But 
when channel blockers are used, we can find out how the 
spiking patterns change. After we focus on finding out the 
relationships between channel blockers and AP change, we 
will examine the mechanism of electrophysiological re- 
sponse to temperature change in neural cells in Aplysia. 

METHODS 

Animals and preparation of the ganglia 

Aplysia juliana and Aplysia kurodai were purchased by 
a local supplier in Jeju-do, South Korea. The species are 
identified according to the literature [16]. 

When an Aplysia was to be dissected, 0.38 M MgCk were 
injected for anesthesia. MgCk, amounting to half weight 
of the Aplysia, were injected. With dissecting tools, the ab- 
dominal ganglia of A. juliana were dissected out and pre- 
served in artificial seawater ASW : L-15 (1 : 1, v/v) as de- 
scribed previously [16]. The abdominal ganglia incubated 
at 34°C for about one hour in 1% protease (type IX, sigma) 



in ASW : L-15. The ganglia treated with protease were 
washed several times with ASW, and were stored at 18°C 
for about five hours. Finally, the ganglion was pinned down 
on a Sylgard plate dish (50x9 mm) and the sheaths of the 
ganglia were dissected out under dissecting microscope to 
expose the neurons on the dorsal side of the abdominal 
ganglion. The medim was exchanged with L-15/ASW. 

Data acquisition 

The Sylgard plate containing the ganglion was placed on 
a small copper plate (50x100 mm) under which two thermo- 
electric controllers were set up (HMN 3940, Acetec Co., 
Korea) side by side. One with the hot side was faced down 
and the other with the cold side was faced down. We in- 
stalled an automatic temperature controller. The data ac- 
quisition (DAQ) card (NI PCI-6221, National Instruments) 
was connected between the Neuroprove Amplifier (Model 
1600, A-M Systems, Carlsborg, WA, USA) and the 
computer. The membrane potential was recorded intra- 
cellularly using the Neuroprove Amplifier with a glass mi- 
croelectrode filled with 3 M KC1. and was sent to channel 
2 of the data acquisition device. For the measurement of 
the temperature near the recording neurons, the PT100 
temperature probe was placed near the Aplysia ganglion. 
To record the temperature, we used a digital thermometer 
(TRM-006, Toho, Japan) connected to a PT100 temperature 
sensor through channel 1 of the DAQ card. The values of 
temperature and membrane potential were simultaneously 
saved onto the computer hard disk. We used the computer 
program previously coded for another experiment. Some 
parts were modified and reused in this experiment. In data 
analysis program, various variables were calculated. Most 
significant variables are the AP amplitude Aap (mV), the 
membrane potential at the negative peak V np (mV), the in- 
terspike interval ISI (ms), the AP half- width duration z)tAP, 
i/2 (ms), the interburst interval IBI (ms), and the sponta- 
neous firing frequency Frequency (s ) [17]. 

Drug application 

Tetraethylammonium chloride (TEA, Sigma: T2265), 
4-aminopyridine (4-AP, TCI: A0414) are purchased from 
Sigma-Aldrich and TCI Co. LTD. By using a pump 
(BT100-2J, LongerPump, Hebei, China), L-15 : ASW (1 : 1), 
10 mM TEA, L-15 : ASW (1 : 1), 10 mM 4-AP, L-15 : ASW 
(1 : 1), 10 mM TEA mixed with 10 mM 4-AP were pumped 
in consecutive order at a rate of 1 ml/min for 30 minutes, 
respectively. This is one cycle. We knew that the concen- 
tration of chemical reagents was at maximum within about 
five minutes before and after each drug treatment, and cells 
were adapted to it. In these time intervals, cells showed 
steady results. So, neurosignals were analyzed using the 
data recorded during these periods. 

RESULTS 

Experimental data for analysis 

To simulate temperature-dependent properties of the 
APs, we performed experiments using Aplysia neurons. In 
these experiments, the cells showed a bursting-oscillation 
pattern at temperatures between 28.8°C and 32.2°C and 
then became silent when the temperature was increased 
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Table 1. Dissection and electrophysiological properties for the experiments A to E and details of the seven different firing patterns 
observed in the acquired data. Silent, beating (bursting with very long burst duration and short interburst intervals), doublet, beating 
chaos, bursting chaos, bursting, and bursting oscillation patterns are shown 



Dissection, electrophysiol. properties and Experimental 



patterns of temperature dependence of APs 


Experiment A 


Experiment B 


Experiment C 


Experiment D 


Experiment E 


Species of animals 


A. Juliana 


A. Juliana 


A. Juliana 


A. kurodai 


A. Juliana 


(weights of animals (g)) 


(220) 


(97) 


(358) 


(316) 


(420) 


Seasons of experiments 


spring 


summer 


spring 


spring 


spring 


Size and positions of cells 


Giant, 


Giant, 


Giant, 


Giant, 


Medium-s., 




d.s. of Icq 


d.s. of Icq 


d.s. of Icq 


d.s. of Icq 


d.s. of Icq 


Media (ext. cell fluid) 


L-15 : ASW 


L-15 : ASW 


L-15 : ASW 


L-15 : ASW 


L-15 : ASW 




(1 : 1, v/v) 


(1 : 1, v/v) 


(1 : 1, v/v) 


(1 : 1, v/v) 


(1 : 1, v/v) 


Range of tem. change (°C) 


10.0-31.3 


11.1-32.4 


13.5-30.1 


12.8-30.5 


6.9-26.6 


Total time for recording (min) 


351 


1,632 


360 


1,079 


2,719 


Total number of spikes 


35,983 


125,487 


35,671 


90,344 


92,257 


Average firing rate (min l ) 


102.5 


76.8 


99.0 


83.7 


33.9 


Sel. dur. for anal, (min) 


35 


55 


29 


66 


59 


Patterns of temp. Silent 


10.9-11.5 


11.1-14.0 


13.5-13.6 




9.6-10.1 


dependence of APs (°C) Beating 


11.5-24.5 


14.0-26.8 


15.0-26.0 


15.2-27.2 


10.1-23.5 


Doublet 


24.5-26.3 


26.8-27.1 


26.0-26.9 


27.2-27.7 


23.7-24.0 


Beating-chaos 


26.3-27.1 


27.1-28.0 


26.9-27.9 


27.7-27.8 


24.0-24.3 


Bursting-chaos 


27.1-28.4 


28.0-28.7 


27.9-28.4 


27.8-28.0 


24.3-24.9 


Square-wave bursting 


28.4-30.0 


28.7-30.2 


28.4-28.8 


28.0-28.9 


24.9-25.9 


Bursting-oscillation 


30.0-31.0 


30.2-32.2 


28.8-30.1 


28.9-31.2 





Ext. cell fluid, external cell fluid; d.s. of Icq, dorsal surface of the left caudal quarter-ganglion; Sel. dur. for anal., selected section 
duration for analysis; Patterns of temp, dependence of APs, Patterns of temperature dependence of APs. 



over about 35°C. However, the temperature had to be re- 
duced to prevent the cell from becoming silent after the 
bursting-oscillation pattern appeared because the silent 
state at high temperature could influence the activity of 
APs. 

Table 1 shows the conditions of each experiment 
(Experiment A to E). In these experiments the species A. 
Juliana and A. kurodai were used. Only large cell somata 
bigger than ~ 100 micron were selected for recording and 
these neurons were located at the dorsal surface of the left 
caudal quarter-ganglion. The temperature was controlled 
between 6.9°C and 32.4°C. The total length of time for re- 
cording was between 351 and 2,719 min. The total number 
of spikes recorded from each neuron was from 35,983 to 
125,487. The average firing rate was between 33.9 min 1 
and 102.5 min \ The duration of recording sections sub- 
jected to data analysis varied from 29 to 66 min. In Table 
1, we selected the five datasets of the experiments that 
showed the following seven spiking patterns sequentially 
as temperature increased: silent, beating (sinusoidal burst- 
ing with a very long burst duration and short interburst 
intervals (IBI)), doublets, beating-chaos, bursting-chaos, 
square-wave bursting, and bursting-oscillation. In these ex- 
periments when the temperature was decreased over a sim- 
ilar range, the above seven patterns reappeared in reverse 
order. For each of the five separate experiments the tem- 
perature intervals corresponding to each spiking pattern 
are presented at the bottom of Table 1, suggesting that the 
firing pattern is likely determined by the instantaneous 
temperature; there were no data acquired below 15.2°C in 
Experiment D (no silent state), and above 25.9°C in 
Experiment E (no bursting-oscillation). 

Fig. 1 displays the averaged values of AP parameters as 
a function of temperature. The IBI has larger errors bars 
than the other AP parameters. The Aap, Atp&, 1/2, and ISI 




^ 10 15 20 25 30 10 15 20 25 30 



T(°C) T (°C) 

Fig. 1. Average values of six AP parameters with error bars. Aap, 
^tAP, 1/2, ISI, and IBI displayed in panels A, C, D, and E, respec- 
tively, decreased as temperature increased. | V np I and Frequency 
displayed in panels B and F, respectively, increased as temperature 
increased. 



decreased with increasing temperature. The IBI in Fig. 1 
E also decreased with increasing temperature; however, the 
error bars in Fig. 1 E from 20°C to 24°C are relatively large 
because different values of the IBI in the sinusoidal burst- 
ing pattern appeared at different temperatures within 
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2000 4000 
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Fig. 2. Temperature and ID as a function of time. Data acquired 
from one temperature cycle in Experiment A are plotted. The near 
symmetry of the graph is taken to mean that the time-dependence 
of the data to be analyzed is reproducible. The small upper box 
displays a magnified view from 3,500 s to 4,000 s. The symbols 
E1A to E4A (Fig. 5), written in the large bottom box correspond 
to the time series displayed in each panel in Fig. 4E1A to Fig. 
5E4A. Sinusoidal bursting with very long burst duration and short 
IBI is shown by the sign looked like "A" in the bottom trace. 



these temperature ranges. |V np | in Fig. IB and the Fre- 
quency in Fig. 1 F increased with increasing temperature. 
Likewise, the spiking frequencies of warm (cold) receptors 
increase as the temperature increase (decrease). 

Based on these experimental results, we wanted to com- 
pose equations to simulate all of the seven different spiking 
patterns. In order to do this, we selected experimental data 
from two of the five experiments, i.e., Experiments A and 
B. Although the selected data had similar firing patterns, 
data from Experiment A displayed more stable firing pat- 
terns, especially in the doublet, beating-chaos, bursting- 
chaos, and square-wave bursting patterns. For this reason 
we chose to use the data from Experiment A to draw multi- 
ple graphs corresponding to each spiking pattern. It was 
also necessary to select a section of the data containing the 
properties of all of the spiking patterns; for this the 
80-minute data which correspond to one temperature cycle 
was selected from the data of Experiment A. The temper- 
ature values and interval duration (ID) of APs are pre- 
sented as functions of time in Fig. 2. In this figure, it is 
shown that when the temperature was increased or de- 
creased over the same temperature range, the same spiking 
patterns reappeared in reverse order. For example, when 
the temperature was decreased from 31°C to 11°C, the sev- 
en firing patterns of Aplysia neurons changed sequentially: 
bursting-oscillation, square-wave bursting, bursting-chaos, 
beating-chaos, doublets, beating, and silent. These spike 
trains in continuous time series are shown in Fig. 3. The 
small box at the top of Fig. 2 displays a magnified view 
of the IDs from 3,500 s to 4,000 s and clearly shows doublets 
appearing. The symbols "E1A"~"E4A (Fig. 5)" written on 
the lower big plot correspond to the names of the panels 
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Fig. 3. The upper and lower traces 
represent temperature and mem- 
brane potential, respectively. When 
the temperature was decreased from 
31°C to 11°C, the seven firing pa- 
tterns of Aplysia neurons changed 
sequentially: bursting-oscillation, squ- 
are-wave bursting, bursting-chaos, 
beating-chaos, doublets, beating, and 
silent. 
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Fig. 4. Comparison of experimental results with results from computer simulations. Panels E1A~E4A display data from Experiment A. 
Panels E1B~E4B display data from Experiment B. The relevant temperatures are shown at the top of each box. Panels S1~S4 display 
values calculated from equations (1) to (4) using parameters from the type 1 model at the temperatures shown in panels E1A~E4A, 
respectively. The variables g K and t u used in the simulations are shown at the top of the boxes in the bottom panel of each set of three 
panels. Panels El A, E1B, and SI display silent pattern. Panels E2A, E3A, E2B, E3B, S2 and S3 display beating pattern. Panels E4A, 
E4B, and S4 display doublet pattern. 
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Fig. 5. Comparison of experimental results with results from computer simulations. Panels E1A~E4A display data from Experiment A. 
Panels E1B~E4B display data from Experiment B. The relevant temperatures are shown at the top of each box. Panels S1~S4 display 
values calculated from equations (1) to (4) using parameters from the type 1 model at the temperatures shown in panels E1A~E4A, 
respectively. The variables g K and r n used in the simulations are shown at the top of the boxes in the bottom panel of each set of three 
panels. Panels E1A, E1B, and SI display beating chaos. Panels E2A, E2B, and S2 display bursting chaos. Panels E3A, E3B, and S3 display 
square-wave bursting. Panels E4A, E4B, and S4 display bursting-oscillation. 
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that are displayed in more detail in Fig. 4E1A to Fig. 5E4A, 
respectively. Several sinusoidal burstings with very long 
burst duration and short IBI are shown by the sign looked 
like "A" in the bottom trace in Fig. 2. Since these signals 
appeared irregularly, we did not simulate these pheno- 
mena. 

The model 

Although the Huber-Braun cold receptor model [14,15] 
describes cold receptors well, proper biophysical model for 
warm receptors has not been searched. To set up nonlinear 
differential equations that could simulate all of these vari- 
ous spiking patterns together, it was necessary to modify 
the equations (Chay-Lee model) derived by Chay and Lee 
[6,7] with temperature-dependent scaling factors; we adapt- 
ed the temperature-dependent scaling factors for the ionic 
currents, p (T), and for the ionic kinetics, <J>(T), that are 
used in the Huber-Braun cold receptor model. To consider 
a temperature-dependent scaling factor for reversal poten- 
tials (or Nernst potentials), a new temperature-dependent 
scaling factor, a (T), was defined and used. The reference 
temperature To was also present in these new equations. 
As a result, we obtained equations using 23 parameters 
with temperature-like scaling factors to simulate the se- 
quential firing patterns shown in Fig. 3. 

The cell membrane of this model contains sodium, cal- 
cium, and potassium channels carrying a fast sodium cur- 
rent, Inh, fast calcium current, Ic a , and a potassium current, 
Ik, respectively. The sodium channel contains a fast-acti- 
vated m gate and a fast-inactivated h gate upon depolari- 
zation. The calcium channel contains a fast-activated d gate 
and a slow-inactivated / gate. The potassium channel is con- 
trolled by an activated n gate. This membrane also contains 
a leak current, II. And then the membrane potential, V, 
is given by 



-C m ^=I Na +I Ca +I K +lL 



= MT)i Na F^t{V-a(T)V N ,} 



1 + e 1 



+ p(T)g c 



-f{V-,7(T)V ca } 



(1) 



1 + e 1 



+ p(T)g K n{V-a(T)V K } 

+ p(T)i L {V-a(T)V L }, 

where C m is the membrane capacitance. The temper- 
ature-dependent scaling factor, p (T), is defined as p(T) = 

T-T„ 

1.3 10 c ; the temperature coefficient of the single-channel 
conductance is Qio=1.35 [18]. The other temperature-depen- 

TH - 2T3 15 K 

dent scaling factor, a (T), is defined as a (T) = 7= — nr7 - - - T ^ . 

1 0 + It a. 15 K 

In Eq. (1), Vnb, Vca, Vk, and Vl are the Nernst potentials 
for the Na , Ca , K , and leak currents, respectively; g Na , 
g Ca , g K , and g L are the maximal conductance for the re- 
spective currents. V m and S m , Vd and Sa are the half-max- 
imal potentials and the slopes at these half-maximal poten- 
tials of the activation variables m and d, respectively. The 
inactivation variables h and f, and the activation variable 
n are all functions of voltage and time and these are de- 
scribed by the following differential equations, 



1 



(2) 



e \ 2S„ I +e \ 2S„ 



1 



£=,(T)^L,( T )^ 
at 7> 



v ; -v \ 
s, I 



(3) 



/v,-v| 


/v-v,| 


J 2S > /. 





Table 2. Seven arrangements of temperature-dependent scaling 


parameters for maximum conductance 


and time 


constants 


Models 


Scaling parameters for 


maximum 


conductances 


Scaling para, for 


Scaling para, for 


Ref. temp. 


Variables 




Na + -current 


K + -current 


Ca + -current CI -current 


the time constants 


Nernst potentials 


CO 




Type 1 


P (T) 


P(T) 


P(T) 


1 


<t(T) 


1 


21 


g K , and t„ 


Type 2 


P (T) 


P (T) 


P(T) 


1 


<t(T) 


1 


21 


go and 


Type 3 


P(T) 


P(T) 


P(T) 


1 


3>(T) 


1 


21 


go and t„ 


Type 4 


P (T) 


P (T) 


P(T) 


1 


<5(T) 


<T(T) 


21 




Type 5 


P(T) 


P (T) 


P(T) 


P(T) 


<S(T) 


1 


23.2 


lo and ~ T h 


Type 6 


P (T) 


P(T) 


P(T) 


P(T) 


<5(T) 


<T(T) 


23.2 


SO and T b 


Type 7 


p(T) 


P(T) 


P(T) 


1 


<t(T) 


1 


11 


go and T l 


Type 8 


1 


1 


1 


1 


1 


1 




go and r h 



Scaling para, for the time constants, scaling parameters for the time constants; Scaling para, for Nernst potentials, scaling parameters 
for Nernst potentials; Ref. temp., reference temperature. 
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dt T i r„ 



(4) 



l + e l 



2S„ 



where the temperature-dependent scaling factor, 4>(T), is 



defined as (J)(T) = 3 10 ° ; many enzyme reactions have a 

Qio value near 3, as does the gating of many ion channels 

[18]. The sigmoidal steady-state inactivations, hoo and fa, 

/ V, - V \ _ 
are defined as hco = l/(l + exp — - ) and foo = l/(l + 

I V<— V \ 

exp — g ), and the steady-state activation, noo, is defined 

/ V - V \ 

as nco=l/(l+exp — ^ ). The parameters h and n ra are 

functions of voltage and time. The relaxation time constants 

V, - V 

r h , t { , and r n are respectively defined as r h = T h /(exp( — - - — ) 

2b h 

V — Vi _ V f -V V — Vt 
+ exp( 2g ' )), r f =r f /(exp(— ) + exp( — )), and r n = 



2S f 



2S f 



V - V 

r„/(l + exp( 2g " )), where r h , 



T f , and r„ are maximal time 



constants; these are also functions of voltage and time. The 
actual forms of these variables are described in equations 
(2) ~ (4). The variables Vh and Sh, Vt and St, V n and S n are 
the half-maximal potentials and the slopes at these 
half-maximal potentials of h, f, and n, respectively. 

We modified the equations derived by Chay and Lee [6,7] 
for excitable cells by adopting temperature-like scaling fac- 
tors defined in the Huber-Braun model. This modified 
Chay-Lee model is required to explain the temperature de- 
pendence of the sequence of the seven different firing 
patterns. Various combinations of p (T), <J>(T), and a (T) are 
shown in Table 2. Scaling factors in the type 1, 2 and 3 
models have similar temperature-dependent scaling factors 
as those in the above mentioned paper by Braun et al. We 
tried to apply scaling factors that depend on temperature 
in this new equations as much as possible. To do this, we 
considered all eight combinations of scaling factors. Table 
2 and 3 contains details of the eight different type models 
obtained by changing the temperature-dependent scaling 

factors. In type 1 model g K and r n were used as variables, 
and the temperature- dependent scaling factors used in the 
Huber-Braun cold receptor model were chosen. Although 
the reference temperature was 25°C in the Huber-Braun 



Table 3. Numerical values of the parameters in the models. The first seven parameters are the same in all eight model types. 
For the middle 11 parameters, numerical values in bold type indicate the parameters that have different values from that used 
in type 1. The last five parameters denote variables 



Numerical values of the parameters in the models 



the model 


LJI11L 




Type 1 


Type 2 


Type 3 


Type 4 


Type 5 


Type 6 


Type 7 


Type 8 


Cm 


u F cm 


-2 


1 


1 


1 


1 


1 


1 


1 


1 




u S cm 


-2 


18 


18 


18 


18 


18 


18 


18 


18 


V Na 


mV 




80 


80 


80 


80 


80 


80 


80 


80 




mV 




140 


140 


140 


140 


140 


140 


140 


140 


V L 


mV 




-71 


-71 


-71 


-71 


-71 


-71 


-71 


-71 


V t 


mV 




-54 


-54 


-54 


-54 


-54 


-54 


-54 


-54 


S m 


mV 




6.4 


6.4 


6.4 


6.4 


6.4 


6.4 


6.4 


6.4 


SNa 


li S cm 


-2 


1,000 


1,000 


1,000 


1,000 


1,000 


1,000 


1,300 


1,000 


V K 


mV 




-90 


-80 


-80 


-80 


-80 


-80 


-80 


-80 


V m 


mV 




-12 


-12 


-12 


-12 


-12 


-12 


10 


-12 


V h 


mV 




-35 


-39 


-39 


-39 


-40 


-39 


-40 


-39 


V d 


mV 




-48.4 


-48.4 


-48.4 


-49.4 


-49.5 


-49.5 


-47.6 


-49.0 


v n 


mV 




15 


15 


15 


15 


15 


15 


14.5 


15 


Sh 


mV 




-5.4 


-5.4 


-5.4 


5.5 


-5.1 


-5.5 


5.1 


-5.5 


Sa 


mV 




5.7 


5.7 


5.7 


5.88 


5.5 


5.7 


5.27 


6 


s t 


mV 




-8.5 


-8.5 


-8.5 


-8.5 


10 


10 


10 


-10 


s» 


mV 




15 


15 


15 


15 


15 


15 


16 


15 


To 


°C 




21 


21 


21 


21 


23.2 


23.2 


11 




So, 


li S cm 


-2 


28 


Variable 


Variable 


Variable 


Variable 


Variable 


Variable 


Variable 










(13-32) 


(15-120) 


(9-20) 


(13-27) 


(7-27) 


(13.2-21) 


(9-27) 


g K 


li S cm 


-2 


Variable 


330 


330 


330 


330 


330 


330 


330 








(120-1,100) 














T h 


ms 




160 


Variable 


200 


120 


Variable 


Variable 


1200 


Variable 










(59-270) 






(123-440) 


(155-230) 




(40-186) 


T n 


ms 




Variable 


80 


Variable 


80 


80 


80 


80 


80 








(40-150) 




(65-700) 












r t 


s 




100 


40 


100 


110 


90 


200 


Variable 


60 




















(12.3-90) 





Para, in the model, parameters in the model. 



378 



NG Hyun, et al 



model, in this model we took a value of 21°C as the refer- 
ence temperature because this represents the middle of the 
temperature range (from 11°C to 31°C) of Experiment A. 
The type 2 and 3 models used the same combination of fac- 
tors as type 1 model except that the variable g Ca , r h was 
used in the type 2 and g Ca in the type 3 instead of g K , 
r n in type 1. The type 4 model is different from the type 
1 because a scaling factor, a (T) is considered in addition 

to the type 1 parameters. g Ca is the only variable used in 
the type 4 model. The type 5 model differs from the type 
1 model because a scaling factor for the maximum con- 
ductance of leak current is considered, and a new reference 



0.2 A 



= 0.1 



0.3 r 



.0.2 



'0.1 



T„/f»(T) 

T n (21*C) = 0.0626 s 



7 1.5 
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o ! 
CO 

E 0.5 
i cb 0 



10 15 20 25 30 



10 15 20 25 30 



10 15 20 25 30 
T(°C) 



CO 1 

E 

£0.5 

P 



10 



-g K (21°C) = 0.536mS CI 



15 20 25 30 

T(°C) 



Fig. 6. Simulation variables. Panel A display the variation in t u 
with respect to temperature. Panel B displays the variation in g K 
with respect to temperature. Panel C displays r n /$(T). In this panel 
the dotted line displays the value of r n /c])(T) when r n is a constant 
value at reference temperature 21°C. The solid line in panel C 
displays r n /$(T) with r n as a variable represented in panel A at 
the corresponding temperatures. At low (high) temperatures the 
value of T n /cJ)(T) is smaller (larger) than the expected value, which 
is shown as a dotted line. Panel D displays p (T) • g K . In this panel 
the dotted line displays the value of p (T) • g K when g K is a 
constant value at reference temperature of 21°C. The solid line in 
this panel displays the value of p (T) • g K with g K a variable 
represented in panel B at the corresponding temperatures. 
Although the dotted line increase as temperature increase, the solid 
line decrease. 



temperature, 23.2°C, is used. The type 6 model is the same 
as the type 5 except for an additional scaling factor, a (T). 
The type 7 model is the same as the type 1 model except 
that the reference temperature is taken to be 11°C and the 
variables g Ca and r f are used instead of g K and r n . Twenty 
one fixed parameters are used in the type 8 model and g Ca 
and r h are variables [19]. 

From these eight types of model, it was necessary to se- 
lect the best type of model to properly explain the temper- 
ature dependence of the sequence of firing patterns. Each 
value of the first seven parameters in Table 3 was the same 
for all models. In the middle section of the table, numerical 
values in bold type indicate the parameter values that dif- 
fer from the corresponding values in the type 1 model. The 
five parameters at the bottom of Table 3 are variables: g Ca , 

Sk> r h> r n > an d r f Sea i s a variable in all of these seven 
types of model except the type 1 model. 

It is, however, necessary to check other merit points of 
the type 1 model. The first merit point is that the type 1 
model was obtained using the same kind of temper- 
ature-dependent scaling factors used by the Huber-Braun 
cold receptor model. The second one is that the temperature 
range of the doublet in the experiment (between 24.5°C and 
26.3°C) was similar to that in simulation results (between 
24.0°C and 26.5°C). The last point to consider is that g K 
and r n were used as variables. These two variables are as- 
sociated with the potassium channels. Interestingly, we 
found that the spiking pattern changed from beating to 
doublet (from Fig. 8E1 to E2), when potassium channel 
blockers, TEA and 4-AP were used. It is noteworthy that 
this change phenocopies the temperature-induced change 
from beating to doublets (from Fig. 4E3A to 4E4A). There- 
fore, the type 1 model was considered to be a suitable simu- 
lation model as it contains the potassium channel asso- 
ciated variables, g K and r n . The procedure to draw figures 
using parameters of the type 1 model is as follows. 

Computer simulation 

The time traces of APs in Fig. 4S1 are plotted from the 
equations (1) to (4) at temperatures indicated in Fig. 4E1A 
using the values of the maximal potassium conductance 
and maximal potassium time constant in Fig. 4S1. Panels 
from Fig. 4S2 to Fig. 5S4 are also composed of the graphs 
plotted from these equations at temperatures indicated in 
the corresponding panels (from Fig. 4E2A to Fig. 5E4A) us- 
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Fig. 7. Ionic currents and para- 
meters as a function of time. Panel 
A, B, and C display Ica(t), lNa(t), and 
lK<t), respectively. The solid, dotted, 
dashed, dash-dotted, dashed-double- 
dotted lines in panels from A to F 
are used to represent the variables 
as a function of time at tempera- 
tures 16.1°C, 21.7°C, 25.4°C, 27.6°C, 
and 31°C, respectively. Panel D 
displays gk(t). The potassium cha- 
nnel relaxation time constants r„(t) 
and the probabilities of openings for 
activation of the potassium channels 
n(t) as a function of time are shown 
in panels E and F, respectively. 
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ing the values of the corresponding potassium maximal con- 
ductances and maximal time constants shown in the panels 
from Fig. 4S2 to Fig. 5S4. In order to draw these graphs, 
we used the following numerical values in the type 1 model: 

(1) membrane capacitance: C m =l V Fcm 2 ; 

(2) the Nernst potentials for ions: VNa=80 mV, Vca=140 
mV, V K =-90 mV, V L =-71 mV; 

(3) the maximal conductance: g Xa =l,000^S cm 2 , g Ca = 
28 vS cm" 2 , invariable (120- 1,100 ^S cm -2 ), i L 
=18 cm 2 ; 

(4) maximal time constants: r h =0.16 s, r f =100 s, r n = 
variable (0.040-0.150 s); 

(5) the half-maximal potentials: V m =— 12 mV, Vh= — 35 
mV, V d =-48.4 mV, Vt = -54 mV, V„=15 mV; 

(6) the slopes at the half-maximal potentials: S m =6.4 mV, 
S h =-5.4 mV, S d =5.7 mV, S f = -8.5 mV, S n =15 mV; 

(7) reference temperature: To=21°C. 

The values of the variables g K and r n used in each simu- 
lation are shown in the bottom panels of Fig. 4 and 5. 

To compare the simulation results with experimental re- 
sults, eight sets of three figures are drawn in Fig. 4 and 
5. The time traces of membrane potentials displayed in the 
top panels of Fig. 4 and 5 (from Fig. 4E1A to Fig. 5E4A) 
show a part of data from Experiment A shown in Fig. 3, 
whereas those in the middle panels (from Fig. 4E1B to Fig. 
5E4B) show a part of data from Experiment B. These sets 
of figures are arranged from low to high temperatures. 
Simulation results of the type 1 model are displayed in the 
bottom panels in Fig. 4 and 5 (from Fig. 4S1 to Fig. 5S4); 
these figures were obtained by solving equations in modi- 
fied Chay-Lee model using Mathematica 5.1. To summa- 
rize, the bottom panels of Fig. 4 and 5 show the simulation 
results of this model at the temperatures shown in top pan- 
els; the used values of variables are shown in bottom 
panels. 

In Fig. 4 panels E1A, E1B, and SI show silent patterns. 
Panels E2A , E3A, E2B, E3B, S2, and S3 show beating, 
whereas panels E4A, E4B, and S4 show doublet patterns. 
In Fig. 5, panels E1A, E1B and SI show beating-chaos, 
whereas panels E2A, E2B, and S2 show bursting-chaos. 
Panels E3A, E3B, and S3 show square-wave bursting, 
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Fig. 8. Comparison of beating and bursting signals influenced or 
uninfluenced by drugs with simulation results. Panel El and E2 
represents time series of APs without drug and with the treatment 
of drugs mixed with 10 mM TEA and 10 mM 4-AP, respectively. 
Panel SI and S2 represents the results of computer simulation with 
modified Chay-Lee model with temperature-like scaling factors. 



whereas panels E4A, E4B, and S4 show bursting-oscilla- 
tions. Although the AP amplitudes in the simulation results 
(panels from Fig. 4S1 to Fig. 5S4) and the corresponding 
experimental results (from Fig. 4E1A to Fig. 5E4A) are 
slightly different, the frequencies of the time series of simu- 
lated APs are all very similar to the corresponding ex- 
perimental results. This is in line with the fact that specific 
intensity information of environmental stimulus is trans- 
formed into a corresponding AP frequency in sensory 
receptors. Therefore, the simulation results of these eight 
model types can be interpreted as similar results to the ex- 
perimental ones despite a little discrepancy in the AP am- 
plitude simulation. 

Fig. 6A and Fig. 6B displays the values of variables, r n 
and g K , respectively, with respect to temperature. Fig. 6C 
displays the maximum time constant r n , divided by the 
temperature-like scaling factor 4>(T), r n /$(T) as a function 
of temperature. The solid line in this figure displays r n /4>(T) 
where r n varies with temperature with each value in the 
panel A; at low temperatures, the values of t u /(J>(T) are low- 
er, and at high temperatures those are higher than the ex- 
pected values, which is shown as a dotted line in this panel. 

Fig. 6D displays the maximum potassium conductance, g K , 
multiplied by the temperature-like scaling factor p (T), i.e., 

P (T) • g K . The solid line in this panel displays p(T) • g K 
where g K varies with temperature with each value in the 
panel B; the values of p(T) • g K decrease, yet the expected 
values of these increase as the temperature increases. 

Fig. 7A, B, C displays the calcium current, Ica(t), the so- 
dium current, iNa(t), and the potassium current, lK(t), as 
functions of time, respectively. The solid, dotted, dashed, 
dash-dotted, dashed-double dotted lines in panels from A 
to F are used to represent the variables as a function of 
time at temperatures 16.1°C, 21.7°C, 25.4°C, 27.6°C, and 
31°C, respectively. As shown in Fig. 7A the maximum cal- 
cium currents increase as the temperature rises until 
21.7°C and then decrease, but sodium currents increase as 
the temperature increases as shown in Fig. 7B and inset 
figure. The maximum potassium current decreases as the 
temperature increases. This is shown in Fig. 7C. Time ser- 
ies of the potassium conductance, gx(t) are shown in Fig. 
7D. The maximum potassium conductance also decreases 
as the temperature increases. Time series of the potassium 
time constants r„(t) are shown in Fig. 7E. The maximum 
values of r n (t) decrease as the temperature rises until 
25.4°C and then decrease as the temperature increases 
further. Fig. 7F displays a time series of the probability 
of the opening of the potassium channel activation gate, 
n(t). In this panel, the maximum values of n(t) at high tem- 
peratures are larger than those at low temperatures. 
However, we need more experimental data to investigate 
the exact mechanisms inherent in determining the temper- 
ature-dependent properties of Aplysia neurons. We then 
carried out the experiments with potassium channel block- 
ers since potassium channels are known to critically regu- 
late the shape of AP and the firing pattern. 

We used two distinct, specific potassium channel blockers 
TEA and 4-AP to find inter-relationships between TEA, 
4-AP and APs. Our electrophysiological recording data re- 
vealed that TEA and 4-AP produced distinct pattern 
changes. Treatment of 10 mM TEA made AP frequency rise, 
width of APs wider, negative peak value rise, and ampli- 
tude of APs decrease. Treatment of 10 mM 4-AP made 
mostly opposite effects. It made AP frequency fall and 
width of APs narrower, and amplitude of APs increase, but 
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negative peak value is not much affected. Mixing both 10 
mM TEA and 10 mM 4-AP showed similar results to 10 
mM TEA treatment (data not shown). We also found that 
the spiking pattern changed from beating (Fig. 8E1) to dou- 
blet (Fig. 8E2) at 19°C when mixed K channel blockers 
(10 mM TEA and 10 mM 4-AP) were bath applied to Aplysia 
neuron. It is noteworthy that this is similar to the change 
from beating (Fig. 4E3A)to doublets (Fig. 4E4A) when tem- 
perature was changed from 21.7°C to 25.4°C. Fig. 8S1 and 
S2 represents the results of computer simulation with 
modified Chay-Lee model using the values of variables 
shown in each panel at 19°C. 

Ionic currents and parameters analyzed in the experi- 
ments with drugs (Fig. 8E1, E2, SI, and S2) or changing 
temperatures (Fig. 4E3A, E4A, S3, and S4) are shown in 
Fig. 9 as a function of time. Fig. 9A and B display lK(t) 
and gK(t), respectively. The solid, dotted lines in the panels 
from A to D shows the time series of ionic currents and 
parameters without drugs and with drugs (10 mM TEA and 
10 mM 4-AP) at 19°C, respectively. The dashed and 
dash-dotted lines in the panels from A to D show the time 
series of those parameters in the firing state of beating at 
21.7°C and in the spiking state of doublet at 25.4°C, 
respectively. The time series of the potassium channel re- 
laxation time constant, r n (t) and the probabilities of the 
openings for the activation of the potassium channels, n(t) 
are shown in the panels C and D. The maximum value of 
r n (t) increases whereas the maximum value of n(t) de- 
creases as temperature increases. 

DISCUSSION 

Modified Chay-Lee model with temperature- dependent 
scaling factors 

By analyzing the equations in modified Chay-Lee model, 
we attempted to trace the mechanism underlying change 




31.55 31.57 31.55 31.57 

Time (s) Time (s) 

Fig. 9. Ionic currents and parameters as a function of time. Panel 
A and B display Ix(t) and gK(t), respectively. The potassium channel 
relaxation time constants r„(t) and the probabilities of openings 
for activation of the potassium channels n(t) as a function of time 
are shown in panels C and D, respectively. The solid, dotted, 
dashed, and dash-dotted lines in panels from A to D indicate 
L15/ASW medium control without drugs at 19°C, L15/ASW medium 
containing 10 mM TEA and 10 mM 4-AP at 19°C, L15/ASW medium 
at 21.7°C, and L15/ASW medium at 25.4°C, respectively. 



of the various spike patterns as a response to a periodic 
change in temperature in Aplysia neurons. 

It is known that the solutions of the dynamic equations 
developed by Chay and Lee [6,7] simulate six consecutive 
spiking patterns with the increasing function of maximal 

relaxation time constant r h (or r n ). Although these equa- 
tions does not contain any temperature-dependent scaling 
factors, the Huber-Braun cold receptor model [14,15] con- 
tains temperature-like scaling factors; the maximal con- 
ductances of currents were multiplied by the temper- 
ature-dependent scaling factor p (T), and the maximal time 
constants were divided by the temperature-dependent scal- 
ing factor <J>(T). To reproduce by simulation the seven differ- 
ent types of signals during an increase (or decrease) in tem- 
perature, a set of equations that contain temperature-de- 
pendent scaling factors are required; the aim was to find 
neuronal equations capable of simulating at once all the 
various spiking patterns shown in Table 1. As a result, we 
crafted a set of equations by modifying the Chay-Lee model 
to include temperature-dependent scaling factors for the 
ionic currents and for the ionic kinetics used in the 
Huber-Braun cold receptor model; computer simulation us- 
ing the parameter values in the type 1 model was reason- 
ably good. However, the simulated AP amplitudes dis- 
played in Fig. 4Sl~Fig. 5S4 were slightly larger than the 
experimental results displayed in Fig. 4ElA~Fig. 5E4A. 
Nevertheless, the frequencies of simulation and experi- 
mental results were similar when compared at the corre- 
sponding temperatures. This demonstrates that we can 
trace the ionic mechanism of the firing pattern change in 
response to temperature change in Aplysia neurons by solv- 
ing these equations of modified Chay-Lee model and ana- 
lyzing the simulation results. 

Ionic mechanism of temperature-driven firing pattern 
change 

The modified Chay-Lee model contains a fast Inh and slow 
lea. The characteristics of these inward currents are similar 
to each other because the d gate in Ic a opens as fast as 
the m gate in Inb. r t controls the burst periodicity, whereas 
r h and r n are related to the amplitude and frequencies of 
spikes [6]. As shown in Table 3, r h , r f , or r n are variables 
in the simulation. These parameters with temper- 
ature-dependent scaling factors decrease as temperature 
increases. This means that the h, f, and n gates open faster 
at high temperature than at low temperature. The max- 
imum conductances g Ca , g Na , and g K with the temper- 
ature-dependent scaling factors theoretically increase as 

temperature increases, whereas g K decreases as temper- 
ature increases as shown in Fig. 6B. The calcium currents, 
Ica(t) have maximum value at 21.7°C and decrease at high 
temperatures as shown in Fig. 7A. In Fig. 7B and C, how- 
ever, maximum values of iNa(t) and Ixft) increase and de- 
crease as temperature increases, respectively. As shown in 
Table 3, we used five variable parameters in the simulation: 
g Ca , g K , T h , r n , and r f . Among the five variables, only two 
are used in the type 1 model: g K and r n ; r n is one of the 
three variables of maximal time constants, and g K is one 
of the two variable of maximal conductances. From these 
figures, it may be supposed that the temperature depend- 
ence of sodium and calcium currents is driven by the tem- 
perature-dependent potassium current. Thus, the volt- 
age-dependent potassium current is assumed to play a lead- 
ing role in the temperature-sensitive reactions of neurons 
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in Aplysia abdominal ganglion. However, more experiments 
have to be performed for further analysis with other chan- 
nel blockers or other currents [20]. 

How can we trace the mechanism of change between the 
various spike patterns in Aplysia neurons in response to 
a periodic change in temperature by analyzing the modified 
Chay-Lee model? To address the issue, we can present the 
following procedure. First, it is necessary to set up non- 
linear differential equations with temperature-like scaling 
factors and experiment results are to be simulated with 
these equations. Second, other experiments should be per- 
formed with ion channel blockers and these experimental 
results have to be simulated with revised equations. During 
the repeated procedure, if all of the AP parameters [17] 
are simulated well, then the problems will be resolved. Let's 
give an instance. In a set of experiments, we found that 
the spiking pattern change from beating to doublet (from 
Fig. 8E1 to E2) at room temperatures in the presence of 
K + channel blockers was similar with the changes from 
beating to doublets caused by change in temperature only 
(from Fig. 4E3A to E4A). Therefore, some of the mecha- 
nisms inherent in the temperature dependencies of Aplysia 
neurons can be expressed by analyzing the results of the 
experiments using modified Chay-Lee model. 

TEA or 4-AP reduced potassium current at a concen- 
tration of 10 mM and blocked ion conduction through potas- 
sium channels in a voltage-dependent manner [21]. However, 
n(t) decreased and r„ increased when TEA was applied [22]. 
These results were consistent with the simulation results 
obtained using modified Chay-Lee model and potassium 
current and conductance decreased, but r n increased as 
drugs (TEA/4- AP) were applied (Fig. 8S1, Fig. 8S2). It 
would be interesting to perform more experiments in the 
future using more channel blockers or at different temper- 
atures with different firing patterns to examine how reli- 
ably our modified Chay-Lee model can simulate ex- 
perimental results. 

Concluding remarks 

When the temperature was increased from 11°C to 31°C, 
the firing patterns changed sequentially from silent to con- 
tinuous spiking, doublets, beating-chaos, bursting- chaos, 
square-wave bursting, and bursting-oscillation sequential- 
ly. We wanted to develop a set of equations to describe these 
firing pattern changes and tried to apply temperature-de- 
pendent scaling factors in modified Chay-Lee model; as a 
result, 23 temperature dependent parameters were used in 
these equations. Figures were drawn by solving these equa- 
tions in this model using Mathematica 5.1. Although the 
simulated AP amplitudes were slightly bigger than those 
in the experimental results, the frequencies of the simu- 
lated and experimental results were very similar to each 
other. These results clearly suggest that our model equa- 
tions could trace the ionic mechanism underlying the firing 
pattern changes in response to temperature change in 
Aplysia neurons by analyzing the solutions of modified 
Chay-Lee model. As a result, the voltage-dependent potas- 
sium current seems to play a leading role in the temper- 
ature dependence of neurons in Aplysia abdominal gang- 
lion. Therefore, it will be of great interest to apply our new 
model equations to simulate the temperature dependence 
of the 14 AP parameters that we identified in Aplysia neu- 
rons in the previous study [17]. 
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